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VORTICES IN THE CO-ORBITAL REGION OF AN EMBEDDED PROTOPLANET 
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ABSTRACT 

We present global 2-D inviscid disk simulations with an embedded planet, emphasizing the non-linear 
dynamics in its co-orbital region. We find that the potential vorticity of the flow in this region is not 
conserved due to the presence of two spiral shocks produced by the planet. As the system evolves, the 
potential vorticity profile develops extrema (inflection points) which eventually render the flow unstable. 
Vortices are produced in association with the potential vorticity minima. Born in the separatrix region, 
these vortices experience close-encounters with the planet, consequently exerting strong torques on the 
planet. The existence of these vortices have important implications on understanding the migration rates 
of low mass planets. 

Subject headings: stars: planetary systems: formation — protostellar disks — hydrodynamics 



planet mass, sound speed, disk viscosity, and the torque 
saturation mechanism. One generic behavior he found is 
that the total torque on the planet first shows an oscilla- 
tion consistent with the librating fluid motion in the co- 
orbital region then it levels off to a constant due to the 
disk viscosity. 

Another study is by Balmforth & Korycansky (2001), 
who used matched asymptotic expansions to explore the 
non-linear dynamics of the flows in the co-orbital region. 
They investigated the role of re-arranging the potential 
vorticity profile as a response to the planet's driving and 
its role in the corotation torque saturation. For inviscid 
flows, they also noticed the formation of vortices which 
they attributed to the secondary instabilities in a forced 
critical layer problem (the critical layer in this case is the 
corotation resonance). 

In this Letter, we present results on the non-linear dy- 
namics of the co-orbital region, based on direct, high res- 
olution numerical simulations (§2). Our approach is quite 
similar to that of Masset's except that we have focused 
on the inviscid limit simulations. As a consequence, the 
total torque on the planet exhibits some interesting new 
behavior, instead of leveling off to a constant as in the case 
of viscous disks. Some of the findings are similar to what 
Balmforth & Korycansky have found in the excitation of 
vortices (§3). We provide an interpretation of why these 
vortices are produced in §4 and discuss the implications of 
our results in §5. 

2. INITIAL SET-UP 

We assume that the protoplanetary disk is thin and 
can be described by the inviscid 2-D Euler equations in 
a cylindrical (r, (p) plane with vertically integrated quan- 
tities. Two gravitating bodies, a central star and a pro- 
toplanet, are located at r = and r = I respectively and 
the 2-D disk is modeled between 0.4 < r < 2. The planet 
is on a fixed circular orbit at r = 1 . The self-gravity of the 
disk is not included. 

The mass ratio between the planet and the central star 
is taken to be ^ = Mp/M^ = 10"''. Its Hill (Roche) 
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1. INTRODUCTION 

The discovery of a large number of short-period extra- 
solar planets has generated renewed interests in the study 
of tidally induced migration of protoplanets in disks. The 
basis of migration is the gravitational interaction between 
a gaseous disk and an orbiting perturber, and has been an- 
alyzed well before extrasolar planets were actually discov- 
ered (Goldreich & Tremaine 1979, 1980; Lin & Papaloizou 
1986a, b). It has been suggested that this tidal interac- 
tion induces the embedded protoplanet to migrate under 
two different circumstances. In Type I migration, the pro- 
toplanet mass is small and the disk structure is weakly 
perturbed. In this limit, the disk response can be ana- 
lyzed with a quasi-linear approximation. Linear calcula- 
tions have revealed an intrinsic imbalance of inner and 
outer disk torques (Ward 1997), so the planet could mi- 
grate quite quickly inward. In Type II migration, the pro- 
toplanet mass is sufficiently large for it to open up a gap in 
the disk (Papaloizou & Lin 1984). The protoplanet's or- 
bital evolution is locked to the evolution of the disk, which 
is likely to be on the viscous timescale (Nelson et al. 2000). 
Both types of migration have been basically confirmed by 
the non-linear numerical simulations. 

The high mobility of small mass protoplanets or embryo 
cores of giant planets from Type I migration has presented 
a challenge to the formation of giant planets. The absence 
of deep oceans is contrary to the inference that any terres- 
trial planet in the solar system may be formed well outside 
the snow-line and migrated to their present location. It is 
thus imperative to explore the cause which may have in- 
hibited the Type I migration. One potential contributor 
to this process is the disk response to low mass planets, es- 
pecially in the co-orbital region where only sketchy analyt- 
ical arguments have been developed. There are, however, 
two notable recent exceptions. Masset (2001, 2002) has 
evaluated the co-orbital corotation torque and found that 
they can play a crucial role for the migration rate. Ex- 
tensive numerical simulations of viscous disks have been 
performed to understand the torque dependence on the 
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radius is th — {fi/'SY^^ — 0.032. The disk is assumed 
to be isothermal with a constant temperature through- 
out the simulation region (i.e., it is attached to a thermal 
bath). The isothermal sound speed, scaled by the Ke- 
plerian rotation speed at r = 1 is Cs/v^.r=i = H/r 
where H is the disk scale height. We used a set of three 
sound speeds Cg = 0.04, 0.05, 0.06. In the simulations pre- 
sented here, we choose an initial surface density profile 
with 5](r) oc r~^/^, in order that the ratio of vorticity to 
surface density (potential vorticity) has a flat radial pro- 
file C = (V X 'v)z/^ ~ canst. (Note that Q has a small 
deviation from being a precise constant due to the finite 
pressure gradient which slightly modifies from its Ke- 
plerian value.) With this choice of initial condition, we 
avoid the generation of inflection points due to the rear- 
rangement of C distribution as it is carried by the stream 
lines. We have also made other runs which have an initial 
potential vorticity gradient and we find that the results 
presented here also hold for them. The detailed analysis 
of these results will be presented elsewhere. 

Since th < H for all three diS^erent sound speeds, the 
planet is embedded. The planet's gravitational potential 
is softened by an approximate 3-D treatment. The key ap- 
proach is to take the surface density at each cell spreading 
it out vertically based on a Gaussian density profile with 
a scale height determined by the local temperature. The 
gravitional force between the planet and the disk matter 
is calculated using this 3-D configuration but all the forces 
with the same (r, (j)) positions are summed together. This 
gives a correction factor C(r, (j)) to the gravitational forces, 
which effectively softens the gravitational potential. The 
torque overestimation in a 2-D geometry was studied in 
D'Angelo, Kley, & Henning (2003) which arises only from 
locations near the planet. As a test if our key finding is 
aff'ected by the 3-D treatment, we excluded several Roche 
radii from the torque calculation and show that it is not 
the case. 

Although there are some accumulation of gas near the 
planet, we do not allow any disk gas to be accreted onto 
it. This boundary condition is equivalent to that of suf- 
ficiently low-mass protoplanets with quasi hydrostatic en- 
velopes (Pollack et al. 1996). The planet's potential is 
switched on over 10 orbits in order for the disk to make an 
adiabatic adjustment. (We have used other, longer turn- 
on times but the results do not depend on them.) 

Simulations are carried out using a 2-D code whose ba- 
sic algorithm was based on FARGO by Masset (2000) and 
Li et al. (2001). The co-rotating frame is used so that 
the positions of the central star and the planet are fixed 
at {r,<f)) = (0,0) and (1,0) (acceleration due to frame 
rotation is also included, see Kley (1998)). Runs are 
made using several radial and azimuthal grids to study 
the influence of resolution, and those presented here have 
{rir X njj) = 300 x 1200. Outflow boimdary conditions are 
used at both the inner and outer radial boundaries. The 
simulations typically last several hundred orbits at r = 1. 

One key difference between our simulations and previ- 
ous studies is that our simulations are performed in the 
inviscid limit, i.e., we do not explicitly include a viscos- 
ity term, though numerical viscosity is inevitable and is 
needed to handle for example shocks. In addition, our 
simulations tend to be of higher resolution. With a radial 



n,. = 300 (some runs go up to n,. = 600, = 2400) and 
0.4 < r < 2, the diameter of the Hill sphere of a planet 
with n = 10~^ is resolved by at least ^ 12 cells in each 
direction. Using a nested-grid technique, D'Angelo et al. 
(2002) have also made very high resolution simulations, 
though their simulations include an explicit viscosity. 

3. TORQUE EVOLUTION 

Our simulation results are generally in agreement with 
previous simulations for the planet on a fixed circular or- 
bit (see Masset 2002) in terms of the generic structures 
produced in the disk by the tidal interactions between 
the planet and its surrounding flows. Since the planet 
mass is relatively small, only density dips are produced 
without gaps (at least during a few hundred orbits simu- 
lated here). We define the co-orbital region here as | Ar| = 
|r — 1| < Qvh- The co-orbital region can be separated into 
several parts depending on their streamline behavior (see 
Masset 2002), including the horseshoe (or librating) region 
|Ar| < rH, the separatrix region th < \Ar\ < and 
the streaming region | Ar| > \fVlrB_. The horseshoe region 
can be further divided into real horseshoe streamlines en- 
closing all three libration points (^3,4, 5) and tadpole like 
streamlines enclosing either Z/4 or L5. 

Our key finding is shown in Fig. 1 which displays the 
total torque on the planet from the disk. It has two 
very distinct phases: Phase I: the total torque is neg- 
ative and smooth, modulated by a period that is sim- 
ilar to the libration period near the L4 and L5 points 
(rjib ^ ■\/4/27/i ^ 40P where P is one orbit at r — 1). 
Phase II: the total torque shows very large amplitude and 
fast oscillations with a quasi-period of a few orbits. We see 
such a behavior in all the runs with three different sound 
speeds, though the times when the "phase" transition oc- 
curs are different. They are at t « 40P, 90P, and 21 OP for 
Cs = 0.04, 0.05 and 0.06, respectively. 

4. SHOCKS, POTENTIAL VORTICITY EVOLUTION AND 

INSTABILITIES 

The obvious question is what causes such large ampli- 
tude and fast changes in torques. Closer inspection of the 
co-orbital region at times around 90P (using = 0.05 
as an example) reveals that vortices are generated in the 
separatrix region. This development is shown in Fig. 2, 
which depicts maps of potential vorticity C through the 
times when the system undergoes its "phase" transition. 
One can see that these vortices are anti-cyclones which are 
also blobs of higher densities (not shown here). Being in 
the separatrix region, they experience very close encoun- 
ters with the planet repeatedly, thus exerting strong "im- 
pulses" to the planet. This results in the large amplitude 
and fast variations in the torque evolution. 

To understand this phenomenon further, in Fig. 3, 
we plot the azimuthally averaged radial potential vor- 
ticity profile (C) for = 0.05 at different times t = 
0, 20P, 40P, SOP, 120P. One sees that this profile deviates 
progressively away from the initially flat radial profile, de- 
veloping two prominent structures at |Ar| « 2 — 'drn, 
with a minimum at Ar « ^Arn and a maximum at 
Ar = 3.4ri/. Comparing with Fig. 2, vortices apparently 
emerge from the minima. 
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The existence of the extrema (or inflection points) in 
the potential vorticity profile is usually regarded as a nec- 
essary condition for non-axisymmetric instabilities in ro- 
tating shear flows (Drazin & Reid 1981). The progressive 
steepening of this profile, say, at Ar w 2.4r// until t < 90P 
and the subsequent flattening (see the curve at t = 120P) 
strongly suggest that there is a secondary instability which 
is excited around Ar » 2ArH and vortices are the non- 
linear outcome of such an instability. 

How could the potential vorticity profile develop such 
inflection points/structures from an initially flat and sta- 
ble profile? As can be seen from the evolution equation 
for potential vorticity 

£)(Ci)/-Di= (C^- V)v = in2D, (1) 
where D/Dt is the Lagrangian derivative, that C is con- 
served along the streamlines in 2D (Korycansky & Pa- 
paloizou 1996). Thus, the initial uniform C distribution 
should be prc;scrvc(i in an inviscid disk. This conclusion, 
however, is only true provided that streamlines can be de- 
fined everywhere in a dissipationless flow. 

Disk flows in the co-orbital region around the planet, 
however, do not satisfy this condition (except perhaps for 
flows within |Ar| < rn) because of the two sets of spiral 
shocks produced by the planet. These two sets of shocks 
emanate from close to the planet and cut through the 
whole disk. Flow lines which encounter these shocks are 
broken and the dissipation in shocks breaks the potential 
vorticity conservation. The consequence of the shock is 
equivalent to that of a strong viscous stress in the momen- 
tum equation, the curl of which, does not vanish. 

To establish the role of shocks in breaking the potential 
vorticity conservation, we have performed three runs with 
different sound speeds Cg = 0.04,0.05,0.06. Even though 
they only span a small range, but the effects are quite 
clear. In Fig. 4, we plot the azimuthally averaged poten- 
tial vorticity profiles for these three runs at a time just 
before the emergence of vortices. Vortices first appear at 
the minimum of the "valley" in the profiles which moves 
outwards as Cg increases. This growth pattern is consis- 
tent with the fact that the starting location of the shocks 
is also moving away from the planet when the sound speed 
increases. 

Note that the cause for the emergence and growth of 
inflection points in the azimuthally averaged C distribu- 
tion presented here is different from that due to its re- 
arrangment along the horseshoe orbits (Balmforth and Ko- 
rycansky 2001). In our calculation with a more general 
initial E (and () distribution, we confirm that onset of the 
shearing instability occurs at an earlier epoch. The ap- 
pearance of the subsequently self-excited vortices is quite 
similar to the development of the Kelvin- Helmholtz (KH) 



instability observed at shear interfaces. A detailed study 
on the properties of this secondary instability, such as its 
threshold and growth rate, will be presented in a forth- 
coming paper. 

5. CONCLUSIONS AND DISCUSSIONS 

In this paper, we consider the tidal interaction between 
a protostellar disk and an embedded planet with a modest 
mass. Linear torque calculations suggest that such inter- 
action leads to angular momentum transfer and inward 
migration of the planets. For planets with mass compara- 
ble to that of the Earth, the inferred migration time scale 
(Ward 1997) is much shorter than the observationally de- 
termined depiction time scale of the disk (Haisch, Lada, & 
Lada 2001). This analytic result raises a problem for the 
formation of giant planets through the core accretion pro- 
cess (Pollack et al. 1996). We reinvestigate this process in 
an attempt to resolve this paradox. 

The embedded protoplanet induces a circulation flow 
pattern near its corotation region in the disk. We show 
here that even in disks with initially uniform potential 
vorticities, inflection points may emerge and grow sponta- 
neously as potential vorticity is generated near the spiral 
shock in the vicinity of the protoplanet. These inflection 
points lead to the onset and growth of secondary instabil- 
ities. 

The existence of secondary instabilities and their non- 
linear outcome as vortices indicate that the flows in the 
co-orbital region are perhaps more complicated than lin- 
ear analyses have suggested. These vortices should exert 
strong torques on the planet. 

One limitation of our current study is that the planet 
is being artificially held on a fixed circular orbit. Such 
strong and rapidly varying torques from the vortices on 
the planet might cause the planet to "oscillate" and lose 
its phase coherence with the surrounding flow. This feed- 
back process should have important implications for the 
Type I migration problem. Allowing the planet to mi- 
grate radially must be included in order to understand the 
planet's response to such torques. These issues will be 
addressed in a forthcoming publication. 

We wish to thank Neil Balmforth, Darryl Holm and 
Edison Liang for useful discussions. This research was 
performed under the auspices of the Department of En- 
ergy. It was supported by the Laboratory Directed Re- 
search and Development Program at Los Alamos and by 
LANL/IGPP. We are also supported in part by NASA 
through NAG5-11779, NAG5-9223, and NSf' through 
AST-9987417. 
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Fig. 1. — Time evolution of the total torque exerted on the planet. Phase I (t < 90P): The total torque is negative, smooth, and shows 
an oscillation period of fa 40P which is consistent with the libration period. We used a turn-on time of 10 orbits for the planet mass. Phase II 
{t > 90P): The torque starts to show oscillations with increasingly high amplitudes. During the transition from Phase I to Phase II, vortices 
emerge around Ar 2.4rjj. These vortices experience close-encounters with the planet, resulting in large amplitude variations in the torque. 
The bold curve depicts a moving average of with a window of 4 = 15P orbits. 
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Fig. 2. — Maps of potential vorticity showing the production of vortices. The planet is located at (Ar, cf>) = (0, 0) and the radial distance 
to the planet is given in units of the Hill (Roche) radius ru- Here, Ca = 0.05. (left) At t = SOP the potential vorticity profile is relatively 
smooth, with two "channels" at roughly |Ar| 2 — 3rj^. (middle) At t = 90P the outer channel shows emerging vortices (anti-cyclones), 
(right) At t = lOOP vortices have grown to their full size in both potential vorticity channels. Note that vortices are moving in the separatrix 
region, so they will experience close-encounters with the planet. 
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Fig. 3. — Evolution of the aaimuthally averaged potential vorticity (C) for Cs = 0.05 at different times t = 0, 20P, 40P, SOP, 120P. Take the 
structure at Ar ~ 2 — 3rH as an example. The profile develops dips (Ar = 2.4rif) and peaks (Ar = SArn) progressively from t = — SOP 
until the vortices are produced at t = 90P, after which the minimum flattens (see the bold t = 120P curve). 
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Fig. 4. — Azimuthally averaged potential vorticity (C) shortly before vortices set in for different sound speeds Cs = 0.04, 0.05, 0.06. The 
evolution for lower Cs is much faster and therefore vortices develop early on. Vortices flrst appear at the minima of the proflle. The locations 
of the minima move outwards proportionally as sound speeds increase, so do the starting locations of the spiral shocks. 



